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Abstract 

Standard methods of structural dynamic 
analysis assume that the structural 
characteristics are deterministic. Recognizing 
that these characteristics are actually statistical 
in nature, researchers have recently developed 
a variety of methods that use this information to 
determine probabilities of a desired response 
characteristic, such as natural frequency, without 
using expensive Monte Carlo simulations. One 
of the problems in these methods is correctly 
identifying the statistical properties of primitive 
variables such as geometry, stiffness, and mass. 
This paper presents a method where the 
measured dynamic properties of substructures 
are used instead as the random variables. The 
residual flexibility method of component mode 
synthesis is combined with the probabilistic 
methods to determine the cumulative distribution 
function of the system eigenvalues. A simple 
cantilever beam test problem is presented that 
illustrates the theory. 


Nomenclature 

r.v. random variable 

CDF cumulative distribution function 

PDS probabilistic dynamic synthesis 

CMS component mode synthesis 

FPI Fast Probability Integration 

FORM First Order Reliability Method 

AMV Advanced Mean Value 

g limit state function 

X vector of random variables 

X* design or most probable point (MPP) 

pf probability of failure 

; u i vector of uncorrelated std. normal r.v.’s 
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{x}’ vector of correlated std. normal r.v.'s 

{x} vector of correlated normal r.v’s 

[C] a Correlation matrix of r.v's in substruc. a 

[L] c a Cholesky decomposition of [C] a 

a number of substructure 

p total number of substructures 

k number of kept modes/substructure 

dot's degrees of freedom 

i internal dof s 

b boundary dot's 

N total dof s/substructure 

Introduction 

Structural analysts have always known that 
the parameters of the system being modeled are 
not deterministic due to manufacturing 
tolerances, material deviation, and other factors. 
Until recently, the primary way to deal with this 
knowledge was to use factors of safety, which 
are qualitative and based primarily on 
experience. In an effort to account for these 
variations in the structural parameters in a more 
quantitative fashion, significant research has 
been performed on developing methods to 
actually use the statistical characteristics of the 
input quantities in the analysis to generate an 
output value that is also described statistically. 
Monte Carlo simulations can be performed to 
calculate these probability distributions, but up to 
a million runs are required for accurate results. 
Approximate techniques have therefore been 
developed that require several orders of 
magnitude less calculations than Monte Carlo 
techniques. One such method, the "Fast 
Probability Integration" (FPI) method 1 , has 
recently been implemented into a new 
probabilistic finite element code, NESSUS. 

Numerical analyses of structural vibration 
generally use the finite element method (FEM) 
as the basis for obtaining free and forced 
response characteristics. A frequent problem 
when using FEM is that, for large models 



composed of many substructures, the number of 
degrees of freedom (dof's) is so large that the 
computational costs are prohibitively expensive 
for eigenanalysis. This problem is particularly 
relevent when using probabilistic techniques 
because the eigenanalysis has to be repeated 
many times. The favored solution method in 
industry for the deterministic situation is to apply 
dynamic component mode synthesis methods 
(CMS), which reduce substantially the dof's for 
the system model. 

This paper defines a procedure for 
combining CMS with probabilistic methods to 
obtain the statistical characteristics in an efficient 
manner. These characteristic are summarized 
in the form of the cumulative distrubution 
function (CDF). The procedure makes use of 
statistical distribution information of each 
substructure's dynamical modes and residual 
flexibility, which are available from experimental 
data. This information is synthesized into a 
system model using the residual flexibility 
method of CMS and the statistics of the system 
dynamic characteristics are obtained using FPI. 
An advantage of this method over existing 
probabilistic structural analysis methods is that, 
in many cases, the statistics of the substructure 
dynamic characteristics may be easier to 
determine than those of primitive random 
variables like geometry, material stiffness, or 
density. Final development of the method 
should allow probabilistic methods to be applied 
to much larger models than previously possible, 
such as turbomachinery bladed-discs, which are 
composed of many almost-identical 
substructures whose structural characteristics 
can be described statistically. 

Prnhahilistic Th f^tiral Rackaround 

Research in the field of probabilistic 
structural mechanics has concentrated in two 
areas. The first can be described as 
perturbation methods, and the second as 
reliability methods. The perturbation method, as 
developed by Collins and Thomson 3 , Kiefling , 
and Collins, Kennedy and Hart 5 , is used to 
derive an analytical expression for the mean and 
standard deviation of structural eigenvalues and 
eigenvectors as a function of the derivatives of 
the statistical characteristics of the mass and 
stiffness matrix for each input random variable. 
Hart and Hasselman 6 used component mode 
synthesis to derive analytical expressions for 
the system, or global, eigenvalue matrix as a 
function of the modally reduced substructure 


stiffness and mass matrices. This approach is 
presently being studied by Mahadevan 7 . 

The research described in this paper 
employs the reliability method approach to 
determining the statistical structural response 
characteristics. To briefly review this technique, 
consider the scalar limit state function g(X), 
defined as 

g(X) = Y(X)-y (1) 

where the limit state value y is some desired 
surface in the random variable space and the 
performance function Y(X) is a function of the 
random variables. This formulation divides the 
space into two parts, g < 0 (Y < y) and g > 0 (Y > 
y). The probability that the function Y not 
exceed the desired value y is the probability that 
g<0. For example, if y equals an eigenvalue of 
interest, p(g<0) would be the probability that the 
actual eigenvalue obtained is less than the 
desired value. 

Cornell 8 , using the reliability approach, 
developed what would later be called the First 
Order Reliability Method (FORM) by truncating 
the higher order terms (hot's) from the series: 

g(X) = g Gix) + ZI 7 - (xi-Hi)+hot’s ( 2 ) 
i=i ° Xi 

where |ix is the mean value of X. This resulted 
in the following simple approximations for the 
mean value and standard deviation of g: 

H g - g(Hx) : (3) 

£ ( H ) °*' 2 ■ (4> 

Now assume that both the original X 
distribution and the resulting g distribution are 
normal, and define g(X) < 0 to be the failure 
region. A transformation of g to standard normal 
coordinates can therefore be performed to 
obtain the probability of failure: 

Z= g(X) - ^ 8 - . (5) 

<>g 

The probability distribution function of a 
standard normal variable has a mean of zero 
and a standard deviation of one, and its 
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cumulative distribution function <U(Z) is tabulated 
in statistics textbooks and software packages. 
Since failure is defined as g(X)<0, then if this 
value is substituted into equation (5), the 
equivalent failure region in standard normal 
coordinates is 


Z < -P (6) 

where the parameter p is defined as 



The probability of failure will therefore be the 
area under the probability density function to the 
left of -B, which is equal to the value of the 
cumulative distribution function <I> at that point. 
Since negative values for p are not tabulated, 
the relationship 

p f = <t>(-p) = 1 -<tKP) (8) 

is used instead to calculate this probability. 

Hasofer and Lind 9 refined this method 
further. They introduced an initial reduction of 
each of the primitive normal r.v.'s xj into 

standard normal r.v.'s vj using 
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In terms of the standard normal r.v.'s vj, the 
joint probability function is bell shaped and 
symmetric about the origin. If the limit state is 
defined to be g(V) = g(vi ,V2), then the minimum 
distance from the line 

g(vi V2) = 0 (10) 


to the origin can be shown to be equal to p. The 
point along the line g = 0 that is closest to the 
origin is called the design point X . It is also 
termed the most probable point because it is the 
point along the line g=0 that has maximum 
probability density. 

The reliability method was expanded by 
Rackwitz 10 to multi-dimensional problems for 
which the limit state curve g=0 is an explicit 
function of the r.v's. Wu and W.rsh.ng 1 
developed the Advanced Mean Value (AMV) 
method, a procedure for using the FORM with a 


minimum number of calculations; this is vital for 
non-explicit limit states, such as finite element 
solutions. The limit state is approximated as a 
linear function about the means of the r.v.’s 
(equation 2) and the partial derivatives are 
approximately obtained by numerically 
differentiating the limit state with respect to each 
r v Values of P and X* are obtained for each 
desired limit state by using the "Fast Probability 
Integration(FPI)" method, which is a compilation 
of the improvements to the FORM made by Wu, 
Rackwitz, and others. At this point in the 
procedure, an exact solution for each of the limit 
states (usually a finite element solution) is found 
by plugging in these most probable points. 
These results and their associated p values can 
then be used to create an entire CDF, which is 
shown in Wu's paper to be very accurate in 
comparison to Monte Carlo simulations for 
several examples. Further iterations can be 
performed by expanding the limit state about the 
new design points instead of about the means 
as was performed in equation (2). The FPI and 
AMV methods have been incorporated in 
NESSUS, a probabilistic finite element program 
presently under development by NASA/Lewis 
and Southwest Research Institute. 

pmhahiiigtin Dynamic Synthesis 

The proposed methodology makes use of 
the residual flexibility method of component 
mode synthesis. This method has been 
developed by MacNeil 12 , Craig and Chang , 
and Martinez 14 . The essential idea in CMS is 
that substructure modes are truncated since 
their higher modes will not have a major effect 
on the system modes. The residual flexibility 
method incorporates the effects of the higher 
modes by determining their flexibility. A side 
benefit is that all the elements of the system 
stiffness matrix can be obtained from test and 
that the mass matrix can be closely 
approximated by a unity matrix in the non- 
boundary partition (equal in size to the number 
of kept modes k). Since all the substructure 
information can be obtained from test, 
probabilistic data can be completely 
incorporated into the system matrices to obtain 
the system modes. 

The first step of the probabilistic dynamic 
synthesis (PDS) method developed in this paper 
is to divide the model of a structure into 
substructures a= a,b,...,p. The physical 
displacement vectors of each substructure, 
which have either a subscript i denoting internal 
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dot’s or a subscript b denoting boundary dot’s, 
can be written as 



where dim xj + dim x& = N total dot’s for that 
substructure. 

Each substructure is represented by n 
samples. Each sample is modally tested 
individually in a configuration such that the 
interface locations with other substructures are 
in a free condition. For substructure cl , sample 


there will be some degree of correlation between 
each of the random variables. These correlation 
values range from zero, or no correlation, to +/- 
1, or fully correlated, and can be easily 
calculated from the measured data. The values 
are placed in a correlation matrix [C] a relating 
each element with every other element. For the 
probabilistic analysis, a set of independent 
standard normal random variables {u} a will have 
to be obtained. This can be accomplished by 
making an orthogonal transformation of [C] a with 
its Cholesky Decomposition lower triangular 
matrix [L] c a to uncouple the {x}’ coordinates, 
thereby creating {u} a 1B . This can be 
expressed for substructures a =a,b,...,p as 


i, the test will yield eigenvalues A? and 

eigenvectors ^ addition, the boundary 

partition of the residual flexibility matrix 
[G w ] a ar e obtained from the measured 
boundary drive point frequency response 
functions of the boundary coordinates 15 . For 
use in the PDS method, only the kept (non- 
truncated) eigenvalues, the boundary 
coordinates of the kept eigenvectors, and the 
boundary partion of the residual flexibility matrix 
are needed. These values can be combined into 

a single vector . defined as 


{xr=[L]“{«r. d3) 

The FPI algorithm requires that each 
independent random variable be varied 
individually by some percentage of its standard 
deviation a, which was chosen to be 50 percent 
for this development, while the other r.v.'s are 
kept constant at their mean values. Each of 
these cases is then back-transformed to form a 
corresponding case of the original correlated 
random variables. These are then plugged into 
the model to generate the limit state 
approximation (equation 2) of the response 
value, which is used to obtain a CDF and the 


M“ = 



\{<P b }k 

^ A 1 ^ 


( 12 ) 


design points X*. 

Since the distributions of the r.v.’s are 
standard normal, .5 a will simply equal a value 
of .5 for the r.v. to be varied. The first case is 
therefore 



,P (14) 


where {(p b } J is a vector of the boundary node 

modal displacements for the j'th mode, and V is 
the j’th eigenvalue of substructure a. 

If the entire sample of substructure oc is 
tested, (x}« can therefore be defined as a vector 
composed of elements that are each a normally 
distributed random variable with measured mean 
and standard deviation. Using equation (9), this 
vector is now converted to {x)’ a , a vector of 
standard normally distributed r.v.'s. In addition, 


The next case will consist of the second element 
in (u} a equaling .5 and all the other elements of 
(u) a as well as all the elements of the other 
{u} m 's equaling zero, and so on. 

For each case the (u) for each substructure 
is then transformed to the set of correlated 
standard normal r.v.'s {x}’ using the transpose of 

[L]“ and then into the original r.v.’s (x) using 
equation (9). The new vectors (X), [4*], snd 
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Full Model 


[Gj are pulled out from {x} and placed in 

substructure mass and stiffness matrices 
according to the residual flexibility formulation: 


K a = 


A + 

sym 


<£ J* 


M a 


7 O' 
0 0 ' 


(15) 



The system mass and stiffness matrices are 
now generated by directly coupling the 
substructure mass and stiffness matrices. This 
is accomplished by ordering the "Kept" dot's of 
each substructure sequentially in the system 
matrices and adding the boundary partitions 
together . These matrices can now be used to 
form the system equation of motion, 


[M} 


sysv 


qsa 

/ qsa 

: 

\ • • 

1 + [K]sys | ' 

qsn 



\ X b 1 

1 *b 



= 0 ( 16 ) 


where the q's are generalized coordinates 
related to the kept modes for each substructure. 


The system eigenvalues are then obtained, 
and a single eigenvalue of interest is chosen. 
This eigenvalue corresponds to a single point on 
the random vector response surface defined in 
the FPI method. As each independent random 
variable in the p number of (u) vectors is varied, 
a new response surface point is obtained. This 
surface can be directly input into the FPI code 
and a CDF obtained for the chosen system 
eigenvalue. This will be only a first 
approximation to the CDF, however. The MPPs 
are plugged back into the substructures' mass 
and stiffness matrices, the system is re- 
synthesized, and new, updated eigenvalue 
levels are obtained for each probability level, 
following the AMV method described by Wu. 
These levels are then plotted to show the entire 
CDF. 


Test Case 

Analysis of a spring-mass system (Figure 1) 
using the PDS method has been completed. 
The test system consists of two substructures, a 


Figure 1 Test Case System 


and b, each having four dof's. 5000 samples of 
each substructure were created initially using 
standard Monte Carlo techniques. To achieve 
complete probabilistic generality, each spring in 
the system was assigned a normal distribution 
with a mean of 200 and standard deviation of 10, 
and each mass was assigned a normal 
distribution with a mean of 1.0 and standard 
deviation of 0.5. The Monte Carlo random 
vectors were then used to create the mass and 
stiffness matrices for the substructures (5000 for 
each) and a modal analysis run on the 
substructure samples to obtain their eigenvalues 
{A) a and eigenvectors [<D]“. Three of the four 
modes for each substructure were "kept" for the 
analysis. The boundary partition of the N x N 

residual flexibility matrix [C w ] a was analytically 

calculated directly from the modes that had been 
chosen to be truncated, in this case just the 
highest one, using 


[G]=I 


i «*+ 1 


MM 

A, 


(17) 


where, in this case, k=3 and N=4. The statistics 
on these dynamic characteristics and the 
correlation between them were then calculated. 
These statistics are listed in Table 1. The listed 
quantities comprise the vectors {x} a and (x} b as 
described in equation (12). 

A distribution characterization routine 17 was 
also performed on the distributions to see if they 
could be characterized as normal, which is an 
assumption of the methodology outlined 
(several more steps would have to be performed 
to handle non-normal distributions). 18 Partial 
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Table 1 Statistics of Dynamic Characteristics 

Substructure A 







30.283 

2.5256 

2 

245.90 

20.897 

3 

553.01 

47.274 

Eigenvectors, boundary location only 

l ocation 


■<3 fcTiT* ft TT» ■ U M l ■ | 

i 

.7072 

.0225 

2 

-.7066 

.05786 

3 


.011097 

Residual Flexibility (one boundary point only) 

Mean 


Standard Deviation 

6.3848 E-04 

3.0688E-04 

Substructure B 






0.0 

0.0 

2 

150.31 

13.402 

3 

488.99 

44.077 

Eigenvectors, boundary location only 



ETi!» F=T¥« 1 ([•]!■ 1 


.53517 

.014106 


.75643 

.043935 


-.66730 

.36929 



Mean 

Standard Deviation 1 


7.4945 E-04 2.981 7E-04 


results of the distribution characterization routine 
for one of the r.v.'s are shown in Table 2. The 
"W" statistic is a good ness-of -fit test developed 
by Wirshing and Carlson, where a smaller 
number indicates a closer fit. 19 The results show 
that the data is well represented by both normal 
and lognormal distributions, with the lognormal 
being a little better. The CDF values for each 
distribution, however, indicate that the curves for 
the two distributions predict almost exactly the 
same value, which can be the case for a 
particular set of lognormal parameters. The 
assumption of normality was therefore deemed 
to be accurate. 

At this point the procedure follows the 
outline discussed previously for an actual case, 
which would use modal testing of physical 
samples to generate the dynamic characteristics 


Table 2 Partial Results of Distribution Types 
Routine 

Substructure B, Eigenvalue 2 
w statistics: (all types included) 


Normal: .00948 

Exponential: .30339 

Weibull: .04448 

EVD: .04323 

Lognormal: .00934 


Normal Distrubution Parameters 
sample mean = 150.61 
sample std. dev. = 6.502 

Normal Distribution CDF fit to data 

response value CDF va l U fl 

143.0 -1210 

147.0 -2895 

150.0 .4628 

155.0 -7503 

Lognormal Parameters, base e 
mu =5.014 
sigma = 0.43201 

Lognormal Distribution CDF fit to data 
response value CDF Yalue 

143.0 -1193 

147.0 -2947 

150.0 -4713 

155.0 .7540 

and correlation matrix. For each substructure, a 
matrix composed of cases of (u) vectors (see 
equation 14) was generated and multiplied by 
the correlation Cholesky decomposition matrix 
[L] c to obtain {x}\ the set of correlated standard 
normal r.v.’s. These were then converted to 
their non-standard normal distributions and used 
in the residual flexibility substructure stiffness 
matrix. The substructures were then coupled 
together and a modal analysis was performed on 
the system matrices. The first system 
eigenvalue for each case, which was the 
response value chosen, was then input along 
with its {u} case into the FPI algorithm routine. 


The output of the FPI routine is the "Mean 
Value Solution," an initial estimate of the CDF of 
the response variable, and the MPP’s for the 
specified CDF probability levels. One MPP from 
the output is shown in Table 3. Following the 
AMV procedure, these MPP's were recorrelated 
and converted to the original dynamic r.v. s as 
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Table 3 Sample MPP Output 

Fund. Evalue Response Value « 0.806648E+01 

Probability 

= 0.010000000 

Most Probable Point 

ftiiha r.v. 

value 

Ut 

-0.5067741 E+00 

u2 

-0.1 808949 E+00 

U3 

0.6107269E-01 

u4 

-0.1579644E+01 

u5 

-0.1287386E+00 

u6 

-0.4253657E-01 

u7 

0.2481 191 E+00 

stih b r.v. 

value 

Ul 

-0.1500363E+01 

u2 

0.2808886E+00 

u3 

-0.2627451 E+00 

u4 

-0.6609528E-02 

u5 

-0.2985903E+00 

u6 

-0.5956755E-01 

u7 

0.22721 07E+00 


before, coupled, and a solution obtained for the 
updated fundamental eigenvalue. This value 
and its associated probability level were then 
used to create a new CDF. 

For verification of the PDS method, a Monte 
Carlo analysis was performed on the same non- 
deterministic spring-mass system with the 
system eigenvalue directly obtained from each 
sample. The CDF for this "full" model is 
superimposed on the AMV CDF from the PDS 
method in Figure 2. A very small amount of 
error is indicated graphically. To identify the 
error quantitatively , the amount of vanation of 
the fundamental eigenvalue from its mean value 
at selected probability values for the PDS 
method was compared to the spread for the full 
model. The result shown in Table 4 indicates 
that the deviations from the mean as computed 
by the AMV and MC methods agree to within 
5%. In addition, the mean value of the 
fundamental eigenvalue computed by the AMV 
method is 8.727, which is only .2 % higher than 
that computed using MC, and the AMV standard 
deviation is .272, which is only 3.2 % less than 
the MC standard deviation of 0.281 . 

fipnrli idinn Remarks 

A new methodology has been presented for 
performing analysis of structures composed of 
substructures whose dynamic characteristics 
can be statistically identified. This method uses 
the substructure eigenvalues, eigenvectors, and 



Figure 2 

CDF's for fund, evalue using PDS 
method and M.C. of full model 


residual flexibility as random vectors for 
determining the desired response value by 
combining new probabilistic analysis techniques 
with the residual flexibility method of component 
mode synthesis. Results for a test case show 
the method predicts close to the same 
cumulative distribution function for the system 
fundamental eigenvalue as a non-substructured 
probabilistic Monte Carlo analysis. 

Future work on this method includes 
examining some basic conceptual questions on 
the limitations and applicability of the method. 
Since some of the dynamic characteristics 
necessary for the synthesis may be difficult to 
measure in some circumstances, the formulation 
of a hybrid method combining analysis and test 
will be pursued. Other questions include finding 
cases where the number of random variables 
can be reduced by perhaps only a|lowing 
stiffness or mass to vary, and examining the 
effect of boundary variability for situations like 
the fir tree interface between blades and discs, 
which are neither fixed nor free. In addition, the 
method will be compared to perturbation 
methods to determine the areas of most efficient 
applicability for each one. 
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8.707 
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0.841 

9.002 


8.987 

0.280 

1 .84% 

0.900 

9.080 

0.353 

9.072 

0.365 

3,35% 

0.950 

9.181 

0.454 

9.169 

0.462 

1 .76% 

0.990 

9.371 

0.644 

9.344 

0.637 

-1.17% 


















































Seminar and Workshop on Reliability Methods in 
Mechanical and Structural Design," given by the 
University of Arizona, Tucson, Arizona, Jan. 20 - 


24 1992 

' 1# Mahadevan,S., "Modern Structural 
Reliability Methods," unpublished report for 
NASA, MSFC, October, 1994, pp. 62-68 

19 Wirshing, P. H., and Carlson, J. R., Model 
Identification for Engineering Variables," Journal 
of the Engineering Mechanics Division, ASCE, 
Vol. 103, No. EMI, February 1977, pp. 125-138 


9 



